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An algorithm for the reduction of massive Feynman integrals with any 
number of loops and external momenta to a minimal set of basic integrals 
is proposed. The method is based on the new algorithm for evaluating 
tensor integrals, representation of generalized recurrence relations jjj for 
a given kind of integrals as a linear system of PDEs and the reduction of 
this system to a standard form using algorithms proposed in j2| , || . Basic 
integrals reveal as parametric derivatives of the system in the standard 
form and the number of basic integrals in the minimal set is determined by 
the dimension of the solution space of the system of PDEs. 



1. Introduction 

Mass effects play an important role in confronting experimental data 
obtained at the high-energy colliders, like LEP and SLC, with theoreti- 
cal predictions. Precise determinations of many physical parameters in the 
Standard Model (SM) require the evaluation of mass dependent radiative 
corrections. In the SM, due to the complicated structure of integrals and 
the large number of diagrams (typically thousands), no complete calculation 
even of two-loop self-energy has been carried out. Many different species of 
particles with different masses have to be taken into account and this makes 
the evaluation of multi-loop integrals a rather difficult problem. Existing 
numerical methods for evaluating Feynman diagrams cannot guarantee the 
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required accuracy for the sum of thousands of diagrams and therefore the 
development of analytical or semi-analytical methods are of great impor- 
tance. 

In calculating Feynman diagrams mainly three difficulties arise: ten- 
sor decomposition of integrals, reduction of scalar integrals to several basic 
integrals and the evaluation of these basic integrals. 

For the reduction of scalar integrals to a minimal set of integrals, recur- 
rence relations are the most suitable tool. Recurrence relations originating 
from the method of integration by parts Q were very helpful in various 
sophisticated calculations. Up to now all attempts to extend the standard 
method of integration by parts to multi-loop diagrams with arbitrary exter- 
nal momenta in the case when all the masses are different were unsuccessful. 
The first successful algorithm for the reduction of two-loop massive propa- 
gator type integrals to a minimal set of integrals || was based on generalized 
recurrence relations proposed in ref. [GJ. 

Using the method described in JD], it is easy to write down an enormous 
number of generalized recurrence relations for any kind of integrals. Not all 
recurrence relations are independent and therefore a prescription for deter- 
mining an optimal set of recurrence relations is highly desirable. Another 
closely related and important problem is the determination of the minimal 
set of basic integrals. 

In the present paper we propose an approach which in principle com- 
pletely solves the problems mentioned before. We formulate an algorithm 
for transforming tensor integrals into a combination of scalar ones. We 
show how to reduce systems of generalized recurrence relations to a linear 
system of PDEs. The algebraic analysis of linear systems of PDEs is an 
intensely developing branch of contemporary mathematics and many exist- 
ing algorithms and methods in this field can be used for our purposes. The 
reduction of linear systems of PDEs to an involutive or standard form or 
to a differential Grobner bases solves the problem of finding the optimal 
set of recurrence relations and a minimal set of basic integrals. Involutive 
system of equations can be used for the analytic or numeric evaluation of 
basic integrals or finding different kind of asymptotic expansion. 

The paper is organized as follows. In Sec. 2 a method for the represen- 
tation of tensor integrals and integrals with irreducible numerators in terms 
of scalar ones with shifted space-time dimension d is shortly described. In 
Sec. 3 we describe the method of generalized recurrence relations. In Sec. 
4 we show how to transform the system of generalized recurrence relations 
into a linear system of PDEs. In Sec. 5 we describe the main features of 
the Standard Form algorithm by G. Reid and explain the correspondence 
between the standard form of the linear system of PDEs and our recurrence 
relations. 
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2. Evaluation of multi-loop tensor integrals 

For the evaluation of L loop Feynman graph amplitudes one needs to 
calculate the integrals 



L . N ni n L 

G^Hum}, K» = n j ddk * n p g mj n n w 

i=l j'=l r=l s=l 



where 

^fc,m = ( k 2 _ m 2 +ie )j ' % = X! ^'"^n + fy'mtfm ( 2 ) 

^ ' n=l m=l 

and g m are external momenta, N is the number of lines, E is the number 
of external legs, u and rj are matrices of incidences of the graph with the 
matrix elements being ±1 or 0. 

The traditional way to evaluate tensor integrals consists of the following 
steps: 

• write for the integral the most general tensor as a polynomial in terms 
of monomials in external momenta and the metric tensor g„ v 

• multiply this tensor by appropriate tensor monomials and obtain a 
linear system of equations 

• solve the system of algebraic equations 

One-loop tensor integrals can always be reduced to a combination of scalar 
integrals without scalar products in the numerator. Evaluating multi-loop 
integrals one encounters so-called irreducible numerators, i.e. scalar prod- 
ucts which cannot be expressed in terms of scalar factors in the denominator. 
The solution of this problem for multi-loop integrals as well as for the eval- 
uation of multi-loop tensor integrals was proposed in [0. The main idea 
of this method consists in the representation of tensor integrals in terms of 
scalar ones with shifted space-time dimension d. For tensor integrals with 
any number of loops, internal lines and external momenta, the following 
formula was proposed: 

N ni n L 



i=l j=l r=l s=l 

L , N 

^ 1 ...A„jte},{^},d+)n j /^n^; >m , (3) 



i=i j=i 
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where T^^ Xn is a polynomial type tensor operator and 







dm 2 ' 



d+G {d) = G {d+2) . 



(4) 



The main ingredients of the derivation of the operator T are indepen- 
dent auxiliary vectors cii(i = 1,...,L) and the use of the a- parametric 
representation. The tensor structure of the integrand on the left-hand side 
of (3) can be written as 



ki^n ■ ■ ■ kix„ 



1 











ini +...+n L 



da 



LXnr 



■ exp [i(a\ki + . . . -I-o^/cl)] 



To convert the integral 

L 



L ,. N 

° (d) = n / ddk * n p z ex p + • • • + a Lk L \, 

%=1 j=l 



Oj=0 

(5) 



(6) 



into the a-representation, we first transform all propagators into a paramet- 
ric form 



i " r° 



da a v 1 exp ia(k 2 — m l + ie) 



[k 2 — m 2 + ie) 1 
Using the d-dimensional Gaussian integration formula 

d d kexp \i(Ak 2 + 2{pk)) 

we evaluate the integrals over loop momenta and obtain 

Ld/2 N 



/ \5 


ip 2 


= *(Sj exp 


~~A 



(7) 



(8) 



G ( d) 



iWio "Jo [D(a)] rf / 2 



x exp i 



uyci) 



i=l 



(9) 



where D and Q are polynomials in a, a and qiqj. Differentiating (9) with 
respect to a,j we get by simple identification the operator T: 



T({q},{d},d 



-iQ({{q},a,{0})p ni o 

- TT _ 



jn\+...+n N 

UL Q 



9ai 
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For any particular tensor integral the operator T can be constructed 
by using computer algebra languages. To generate T for 2-3 loop tensor 
integrals of rank 3-4 with the help of a short package written in FORM, 
several minutes on a PC Pentium 90 were needed. To our opinion the 
evaluation of tensor integrals by the proposed method has several advantages 
in comparison with the traditional method. Firstly, in order to obtain the 
tensor decomposition no contractions with external momenta and the metric 
tensor and no solution of a linear system of equations are needed. Secondly, 
it is easy to select the scalar coefficient of the particular tensor structure. 
Thirdly, a representation in terms of integrals with shifted d (as we have 
seen on many one- and two-loop examples) is very compact and may be 
useful for numerical calculations. 

To obtain the final result, scalar integrals with different indices and 
different shifts in d are to be evaluated. This problem can be solved by 
using the method of generalized recurrence relations as proposed in Q. 

3. Generalized recurrence relations 

In order to obtain recurrence relations for integrals we use identities like: 

where R is an arbitrary tensor polynomial. After performing the differenti- 
ation, two different representations for scalar products can be used: 



a) hqj = ~(kf + q] - - qj) 2 ), 




By using all possible combinations of these representations, we produce 
many relations connecting integrals with changed exponents of scalar propa- 
gators and changed values of the space-time dimension. For the same scalar 
product both a) and b) can be used. Combining the different relations, one 
can try to find the most optimal set of relations for the reduction of the 
concrete class of integrals to the minimal set of basic integrals. 

In fact, in the traditional method of integration by parts only represen- 
tation a) for scalar products is used. Our derivation is more general and 
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it includes the integration by parts method j|] as a particular case. Recur- 
rence relations connecting integrals with shifted d cannot be obtained from 
the traditional method of integration by parts. 

To illustrate the difference, we consider the one-loop propagator type 
integral with massive particles: 

T(d) - f ddkl p»i p"2 /io\ 
1 flV2 J ^ 7r d/2]- r fcl,0- r fci-g,m- \ la ) 

From the traditional method of integration by parts two relations can be 
derived: 

2u 2 m 2 I^ V2+1 + u 2+ i + »l(rn 2 -q 2 )li? +1 U2 

-(d-2v 2 -v 1 )lW 2 = 0, 

^"Wl ^-1-^2^-1 u 2 +l + Vl\ m ~Q ) I u 1 +l u 2 

^ 2 ( m 2 + q 2 )ll d \ 2+1 + (u 2 -u 1 )llf U2 = 0. 

(14) 

The integral lj)fl 2 (q 2 , 0, m 2 ) is proportional to the Gauss hypergeometric 
function QJ: 

(m2r+"-4r(|)r(»i) 



^2 

7/, 7/-, -U jy„ _ • , 

X 2 ^l 



z/i,/vi +i/ 2 - f ; (f_ 



(15) 



As is well known there are fifteen relations of Gauss between contiguous 
functions 2F1 ■ Substituting (15) into (14) one can find correspondence be- 
tween the recurrence relations (14) and only six relations of Gauss. The 
reason is obvious - in the relations (14) the third parameter of 2F1 in (15) 
does not change and therefore not all corresponding relations for contigu- 
ous functions can be reproduced. If we include into consideration also the 
generalized recurrence relations, we cover all fifteen relations. 

Usually the number of generalized recurrence relations for a given kind 
of integrals is quite substantial. For their effective applications we must find 
answers to the following questions: 

• How can one use the generalized recurrence relations? 

• How to find the minimal set of recurrence relations? 

• How to find the minimal set of basic integrals or how to determine the 
number of basic integrals? 
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We propose to solve these problems by transforming the system of general- 
ized recurrence relations into a (in general overdetermined) system of linear 
partial differential equations (PDE) and we apply algebraic methods being 
developed by mathematicians for their analysis. Several algorithms and dif- 
ferent computer algebra packages for the analysis of linear systems of PDE 
exist. In our opinion the most adequate will be the algorithm "Standard 
form" by G. Reid 



4. Transformation of generalized recurrence relations into a 

system of linear PDEs 

The system of generalized recurrence relations can be transformed into 
a linear system of PDEs by several methods. Two procedures turn out to 
be most convenient. 

In the first approach recurrence relations for integrals with arbitrary 
powers of propagator indices Uj are considered. If a term with i>i — 1 occurs 
in the recurrence relation then one should make the substitution 

Vi—*Vi + 1- 

Thus we will obtain equivalent recurrence relations connecting integrals with 
positive shifts of indices. Propagators with positive 'shifts' of indices must 
be represented as: 



i T{n) d r 1 



[(fci - q s ) 2 - mfp+ r T{ui + r) (dm 2 ) r [fa - q s ) 2 - m 2 ] 1 



(16) 



Integrals I with different shifts of the parameter of the space-time dimen- 
sion d must be considered as different functions, i.e: 

I d = V\ I d+2 = V 2 , I d+4 = V 3 ,... (17) 

These substitutions allow one to transform the system of recurrence relations 
into a linear system of PDEs for the vector function V = {V 1 ,V 2 , . . .}. 
Such a transformation of recurrence relations into a linear system of PDEs is 
convenient for the analysis of integrals with arbitrary indices. Unfortunatly, 
first approach leads to essential technical difficulties. 

In the second approach one considers a system of recurrence relations 
for integrals with particular integer values of vi > 0. First, one should 
derive recurrence relations for the given integral lj)f\.v N with arbitrary i/'s. 
From relations obtained, one needs to derive recurrence relations for the 
particular sets of integer values of v% > by increasing the sum of indices 
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S = N, N + 1, N + 2, ... (S = J2i=i u i)- m this case the system of linear 
PDEs will include integrals corresponding to graphs with a different number 
of lines. These integrals are to be considered as different functions. As in 
the previous case scalar propagators with Vi > 1 must be represented as 



[(fc _ q s )2 _ m 2]r+l r(r + 1) (am 2)r [(fc. _ ^)2 _ m 2] > V ^ 

and integrals with shifts in d must be considered as different functions. 

As was mentioned before the system of linear PDEs for a given integral 
will include the original as well as simpler integrals obtained by contracting 
lines in the original one. For integrals with contracted lines one should 
also write the system of equations. For all these integrals one can consider 
different sets of i/'s which can be classified according to the value of S. 

The analysis by means of "Standard Form" should be started with the 
system for the simplest non-zero integrals obtained from the original in- 
tegral by contracting as many lines as possible. Applying Standard Form 
algorithms the minimal set of recurrence relations (leading derivatives) and 
minimal set of basic integrals (parametric derivatives) for the simplest inte- 
gral will be found. This information will be used to investigate more com- 
plicated integrals having more lines. Using the terminology of Q integrals 
with contracted lines can be treated as classification functions. 

The second approach to our opinion is more suitable for considering 
integrals occurring in calculating Feynman diagrams. In the first approach 
with arbitrary non-integer indices one can expect more integrals in the basic 
set. 



5. Algorithm Standard Form 

As we already mentioned the algorithm Standard Form (SF) was pro- 
posed by G. Reid in [Q, |3| and was implemented as Maple package in M. In 
this section a short review of the SF algorithm will be given. More detailed 
information can be found in g |, §. 

In order to compare different differential equations and different terms 
in a differential equation one should introduce total ordering which will be 
denoted as > s . The concept of total ordering is very important. 

Let us introduce several shorthands: 

a = (ai, . . . , a m ), 

flai+...+a m 

n vp — yv 

dx?...da%r ' 

ord(a) = ai + . . . a m > 0, order of the derivative (19) 
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The most complete list of different orderings the reader can find in 0]. 
By default in the SF package Tresse ordering is used. We say that 
D a V p > s D b V q if 

(i) ord(a) > ord(b) 
or (ii) ord(a) = ord(b), and p < q, 
or (hi) ord{a) = ord(b), p=q and a >i ex b, 

where the lexicographical ordering >i ex is defined by a >\ ex b iff the leftmost 
nonzero a& — bk > 0. Other orderings can be used in the SF package. 

For the reduction of Feynman diagrams to the minimal set of basic 
integrals we found that the following ordering is more efficient than Tresse 
ordering: 

(i) p > q 
or (ii) p=q and ord(a) > ord(b), 
or (hi) ord(a) = ord(b), p=q and a >i ex b. 

This ordering allows one to exclude all functions with shifts in d from 
the minimal set of differential equations. 

Input for the SF algorithm may be any linear system of PDEs for the 
vector function V = {V , ...,V P } depending on m independent variables 
x%, ...,x m . In our case V q will be a scalar integral, the same integral with 
different shifts in d and integrals obtained from the given integral by con- 
tracting lines. As was mentioned in the previous section integrals with 
contracted lines can be treated as classification functions. As independent 
variables x%, ...,x m one can take masses or their ratios. We assume that all 
lines of the diagram have different masses. 

The main result of the application of the SF algorithm will be a set of 
leading derivatives 

D h V v = fl (20) 

where /? are explicitly given functions of xi and derivatives of V q such that 

i) the derivative on the l.h.s. of (|20| ) is strictly higher in ordering > s 
than any derivative on the r.h.s. of (|20|) 

ii) no derivative appears on both the l.h.s. and r.h.s. of (^C|) 

hi) the derivatives on the l.h.s. and r.h.s. of ( |20| ) are all distinct 

iv) no derivative in ( p0[ ) is a nontrivial derivative of any derivative in 
the l.h.s. of ©. 

v) the integrability conditions of ( |20| ) are identically satisfied modulo all 
lexicographic substitutions which follow from ([20]). 

To clarify the last statement we remind the reader of the definition of 
integrability (consistency) conditions. For any distinct pair of equations 

D a V p = fl D b V p = fl (21) 
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and for any c = (c^) £ N m , > max{afc,6fc} the consistency conditions 
are: 

D c V p - D c V p = D c ^Jl - £> c _ 6 /* = 0. (22) 

Derivatives on the r.h.s. of ( p0|) are called parametric derivatives i.e. by def- 
inition a parametric derivative DhV v is a derivative that cannot be obtained 
by differentiation of any leading derivative. 

The standard form of a linear system of PDEs is achieved by repeating 
the following process until all conditions of the standard form are satisfied: 

i) isolate, and solve for, the highest order (leading) derivatives of each 
equation 

ii) substitute back from (i) throughout the rest of the system 

iii) append new equations resulting from the integrability conditions. 
The output of the Standard Form algorithm (20) can be identified with 

our recurrence relations as follows: 

• the set of equations for leading derivatives corresponds to a minimal 
set of recurrence relations 

• parametric derivatives are our basic integrals 

• the dimension of the solution space of a system is equal to the number 
of different parametric derivatives in the r.h.s. of (p0|). 

To prove finiteness of the dimension of the solution space we refer to 
the known fact that Feynman integrals are holonomic functions. Higher 
order derivatives (integrals corresponding to scalar diagrams with dots on 
lines and/or shifts in d) can be obtained by differentiating the set of lead- 
ing derivatives (p0|). In the Standard Form package Q) there is a special 
subroutine for reducing higher derivatives to parametric ones. Thus, in 
principle, after tensor reduction this subroutine can be directly applied to 
the evaluation of Feynman diagrams. 

Using the Standard Form package we performed classification of the 
generalized recurrence relations for the one-loop propagator type integrals 
and for the simplest two-loop vertex integral with four lines. 

Unfortunately the Standard Form package ||] (written in Maple) is not 
powerful enough and cannot be used for the calculation of real diagrams in 
gauge theories. 

We believe that the implementation of the Standard Form algorithm in 
other computer languages will give better performance and can be used in 
practical calculations. 

In the present paper we gave only a short description of the general 
scheme for evaluating arbitrary Feynman diagrams. Details and develop- 
ments of the algorithms presented will be published elsewhere. 
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